Simulation of a stationary dark soliton in a trapped zero-temperature Bose-Einstein 
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We discuss a computational mechanism for the generation of a stationary dark soliton, or black 
soliton, in a trapped Bose-Einstein condensate using the Gross-Pitaevskii (GP) equation for both 
attractive and repulsive interaction. It is demonstrated that the black soliton with a "notch" in the 
probability density with a zero at the minimum is a stationary eigenstate of the GP equation and can 
be efficiently generated numerically as a nonlinear continuation of the first vibrational excitation 
of the GP equation in both attractive and repulsive cases in one and three dimensions for pure 
harmonic as well as harmonic plus optical-lattice traps. We also demonstrate the stability of this 
scheme under different perturbing forces. 

PACS numbers: 03.75.Lm, 05.45.Yv 
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I. INTRODUCTION 

Solitons are solutions of wave equation where localiza- 
tion is obtained due to a nonlinear interaction and have 
been observed in optics Q , water waves [| , and in Bose- 
Einstein condensates (BEC) 0ildjl3- The bright solitons 
of BEC represent local maxima p and the dark and grey 
solitons local minima [E IS • A stationary dark soli- 
ton where the local minimum goes to zero value is called 
a black soliton. There have been experimental study of 
bright Q, dark and grey [1 Q solitons of BEC. Dark 
solitons of nonlinear optics |lj are governed by the non- 
linear Schrodinger (NLS) equation which is similar to the 
mean-field Gross-Pitaevskii (GP) equation [T(j describ- 
ing a trapped BEC. More recently, dark solitons have 
been observed in trapped BECs 0, Q ■ 

The one-dimensional NLS equation in the repulsive or 
self-defocusing case is usually written as 



iu t + u xx - \u\ 2 u = 0, 



(1) 



where the time (t) and space (x) dependences of the wave 
function u(x, t) are suppressed. This equation sustains 
the following dark and grey solitons 



u(x, t) = r(x — ct) exp[— i{<f>(x — ct) — fit}], (2) 



with 



r 2 (x-ct) = 7 ? -2^ 2 sech 2 [^(a;-ci)], (3) 
<j>{x-cb) = tan" 1 [-2^/ctanh{^(x- ct)}], (4) 
£ = V&V ~ c*)/2, (5) 



where c is the velocity, fi is a parameter, and rj is related 
to intensity. Soliton having a "notch" over a back- 
ground density is grey in general. It is dark if density 



trapped zero-temperature BEC Q . It has been suggested 
that the black soliton of a trapped BEC could be a sta- 
tionary eigenstate of the GP equation 0,0 as in the case 
of the trap-less NLS equation. Here we re-investigate the 
origin of the black soliton in a trapped zero-temperature 
BEC and point out that this soliton 0, H El H 13 
is the first vibrational excitation of the GP equation for 
both attractive and repulsive atomic interactions and is a 
stationary eigenstate. We suggest a scheme for numerical 
simulation of a stationary dark soliton by time evolution 
of the linear GP equation starting with the analytic vi- 
brational excitation, while the nonlinearity is slowly in- 
troduced. We simulate a stationary dark soliton in a 
harmonic and harmonic plus optical-lattice traps in one 
and three dimensions. In all cases the simulation pro- 
ceeds through successive eigenstates of the GP equation. 
Consequently, the stationary dark soliton in a trapped 
BEC could be kept stable during numerical simulation. 
To illustrate the stability of our scheme we also study the 
breathing oscillation of the stationary dark soliton upon 
application of different perturbations. 

The stationary dark soliton being an excited state is 
thermodynamically unstable. It is also unstable due to 
quantum fluctuations However, these instabilities 

are not manifested in the mean-field model. It seems that 
it will be difficult to generate black solitons experimen- 
tally because they are fragile to perturbations, transverse 
dimensions, quantum effects, and thermal perturbations, 
etc. Nevertheless, we demonstrate that they are station- 
ary eigenstates of the GP equation and exploit this infor- 
mation to illustrate a simple numerical scheme for their 
generation. The present numerical scheme has recently 
been used successfully to simulate the stationary dark 
solitons of a degenerate boson-fermion mixture 
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at the minimum. At zero velocity the soliton IL BLACK SOLITON IN A HARMONIC TRAP 



becomes a dark soliton: |it(x, i)| = ^/rj tanh[y/ {rj/2)x]. 

The similarity of the NLS equation |JTJ to the GP equa- 
tion of a trapped BEC (Eq. JJJJ below) imply the possi- 
bility of a stationary dark soliton, or a black soliton, in a 



The mean-field dynamics of a trapped BEC is usually 
described by the time-dependent GP equation 0] . For a 
strong radial confinement in an axially-symmetric config- 



2 



uration, the GP equation can be reduced to the following 
quasi-one-dimensional form [l2l IT3L fl5| 

iu t + u xx - n\u\ 2 u — V(x)u, (6) 

where a positive nonlinearity n represents repulsive (self- 
defocusing) interaction and a negative n represents at- 
tractive (self-focusing) interaction. In Eq. V(x) is 
the external trapping potential. The normalization of the 
wave function is given by /_„ | w | 2 c?ar = 1. The reduction 
of the GP equation from three to one dimension can be 
performed in a straightforward fashion for a single- |l8| 
as well as coupled-channel [Tflj | cases for small nonlinear- 
ity. Nevertheless, for large nonlinearity corrections are 
needed [20|. However, we shall neglect these corrections 
here. 

There is no known analytic solution to Eq. JHJ for 
V(x) 7^ and n ^ 0. Of course, for n — we have the 
well-known harmonic oscillator solution for a harmonic 
trap. However, for V(x) = and positive n, Eq. © has 
the following unnormalizable dark soliton: 

u(x,t) = \/2/ntanh(x) exp (— 2it). (7) 

Soliton J2J) has a stationary notch with zero minimum at 
x = on a constant background extending to x — ±00. 

It has been conjectured that a stationary normaliz- 
able dark soliton exists in Eq. 10 with a harmonic trap: 
V(x)= x 2 and satisfies the same boundary condition as 

HBEI3 

u Ds( x ) — N t&nh(x)us(x), (8) 

where itrjg(x) is the black soliton, ug(x) the ground-state 
solution to Eq. (JSJ) and N the normalization. The form 
(JSJ has been used [13j as an initial guess to a fixed point 
algorithm that finds the exact numerical stationary solu- 
tion. 

The ansatz (JSJ has been used as starting guess in nu- 
merical studies on dark solitons. Actually, in some appli- 
cations |13| the Thomas-Fermi (TF) approximation u^-p 
|10| has been used in place of us for positive n: 

u s (x) w u TF (x) = v/max(0, j/t - V(x)]/n), (9) 

where max(,) denotes the larger of the two arguments 
and fx is the chemical potential for solution us(x). 

The ansatz (|SJ| is not an eigenstate of Eq. JSJ . Assum- 
ing that it is close to an eigenstate, numerical iteration 
of Eq. JfJJl should lead to a stationary dark soliton at 
large times. Whenever the input JHJ is not a good ap- 
proximation to the stationary dark soliton, oscillations 
are expected upon iteration. Kevrekidis et al. fl3l ] give 
a detailed parametric linear stability analysis of the sta- 
tionary solution obtained through their fixed point iter- 
ation scheme varying different parameters, and find win- 
dows of stability and windows of instability. They also 
extend their discussion on instability to the case of an 
optical-lattice potential. In many cases no stable soliton 
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FIG. 1: The dark soliton |u(x)| of Eq. © with V(x) = x 2 
vs. x for n = 1 obtained by time evolution with input © 
at t = 0, 10000, 20000, 30000, and 40000. The dark soliton 
oscillates as a grey soliton without ever converging. 



has been obtained, or an oscillating grey soliton has been 
found 0. 

One way to achieve a stable soliton with initial guess 
©, which is not an eigenstate of Eq. (JBJ, is to include 
a dissipation in the system. The numerical iteration of 
a slightly inaccurate solution would generate radiation 
in general and usually not converge to any stationary 
state without dissipation. This could be the source of 
instability in Ref. [l3| ■ In the following, first we study the 
deficiency of using Eq. (JHJ in Eq. JJjJ without dissipation 
in generating a stationary dark soliton and then illustrate 
the present alternative scheme, where we use an exact 
eigenstate of Eq. © to generate the stationary dark 
soliton. 

First we performed extensive calculations using ansatz 
(JSJ) in the time evolution of Eq. iJBJ using the Crank- 
Nicholson algorithm [2lj for different n. We discretize the 
NLS equation with space step 0.05 and time step 0.0025, 
which was enough for achieving convergence of a station- 
ary problem. We used accurate numerical solution to 
us(x) in place of the TF approximation (jj|J|. The present 
Crank-Nicholson algorithm is appropriate not only for 
the calculation of stationary states but also for nonequi- 
librium dynamics pH with absorptive potential during 
collapse |22j . 

In place of the Crank- Nicholson numerical scheme used 
in this paper one could also use the method of imagi- 
nary time propagation |23| for obtaining the stationary 
states of the GP equation. By replacing the time by an 
imaginary time variable, the original time-dependent GP 
equation becomes a diffusion-like equation, and propa- 
gation in imaginary time leads to relaxation towards the 
vibrational ground state. The imaginary time propagator 
can be expanded in Chebychev polynomials which leads 
to a stable and efficient scheme. This approach can be 
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FIG. 2: The dark soliton |u(a;)| of Eq. @ with V{x) = x 2 
vs. x for n — (a) 100 and (b) —3 obtained by time evolution 
with input (JHJ at different times t. Also plotted are present 
stationary results obtained with input 1101 marked "pr" (full 
red line) which do not change with time. 

easily modified to obtain vibrationally excited states as 
well. By repeating the relaxation (imaginary time prop- 
agation) but filtering out any contribution of the ground 
state at each time step with anti-periodic boundary con- 
ditions, one obtains the first vibrationally excited state. 
This could be an interesting future work. 

Now we perform a direct time evolution of the full GP 
equation with ansatz (jSJ as input and show in Fig. 1 the 
solution at different times for n — 1 . The authors of Ref. 

are making a different time evolution, e. g., a per- 
turbation of their exact stationary solitary wave with a 
uniformly distributed random field of amplitude 0.01. As 
the successive states are not stationary eigenstates these 
schemes may run into numerical difficulty when the non- 
linearity is large. Even for a relatively small nonlinearity 
of ri = 1, in that approach the solution does not con- 
verge at large times: the initial black soliton becomes a 
grey soliton and oscillates around a mean position at the 
center of the trap. Although the wave-function density 
is symmetric at t = 0, it becomes non-symmetric with 
the evolution of time. This makes the black soliton to 
oscillate as a grey soliton upon time evolution before be- 
ing destroyed eventually for much larger values of time 
t. This trouble as noted in Fig. 1 increases with the 
increase of nonlinearity n. 

To circumvent the above-mentioned problem, we find a 
direct solution to Eq. ijfJJl for the stationary dark soliton 
with the asymptotic boundary condition implicit in Eq. 
©, e. g., Mrjg(x) ~ x as x — > and itTjg(x) — * as 
x — > ±oo. The solution satisfying these conditions is the 
nonlinear evolution of the first vibrational excitation of 
the linear oscillator, obtained by setting n = in Eq. 
©: 



ui{x,t) = J (2/V^r)a;exp(-.x 2 /2)exp(-3it). (10) 

The possible stationary dark soliton of Eq. JSJl can be ob- 
tained by time evolution of the GP equation with u\ (x, 0) 
as input at t — 0, setting n — 0. During time evolu- 



0.5 




x 



FIG. 3: The stationary dark soliton |m(i)| of Fig. 2 (a) 
at times t = 0,10000,20000,30000,40000,50000 calculated 
using the present scheme with input H0H as n is suddenly 
changed from 100 to 120 at t = 0. 



tion the nonlinearity n should be slowly introduced until 
the desired nonlinearity is achieved. In this work we in- 
creased n by 0.001 at each time propagation, which was 
sufficient for convergence. By this procedure a stationary 
dark soliton could be obtained for very large n. 

Next we compare the time evolution of the dark soliton 
using conventional ansatz JSJ and the present suggestion 
based on Eq. (|10fl . The results of numerical simulation 
using the two schemes are plotted in Figs. 2 (a) and (b) 
for n = 100 and —3 at different times. We show results 
for these two n, as we found that convergence was more 
difficult for negative n and large positive n values. The 
case n = —3 discussed here is specially interesting as we 
demonstrate, contrary to popular belief, that dark soli- 
tons of a trapped quasi-one-dimensional BEC could be 
stationary for attractive interactions as well. The iter- 
ative solution using Eq. (jSJ) may execute oscillation on 
time evolution, whereas the solution from input l|10l) re- 
mains stationary on time evolution as the system passes 
through successive eigenstates and results in a stationary 
dark soliton of the full NLS equation. In Figs. 2 we show 
the result for the soliton using the present procedure only 
at t = as this result does not change with time. The re- 
sult based on Eq. (JSJ oscillates on time evolution as can 
be found from Figs. 2. However, the oscillation is not so 
severe for small repulsive n (not shown). The oscillation 
increases for large repulsive n as well as for attractive 
nonlinearity. From Figs. 1 and 2 we see that for small n 
the oscillation is severe for t > 10000 whereas for large 
n(= 100) it is disturbing even for t = 100. From Figs. 2 
we find that a direct numerical solution for the first ex- 
cited state of Eq. © is the stationary dark soliton that 
we look for. 

The fact that a certain numerical scheme converges to a 
stationary state may not necessarily signify its stability. 
Although the dark solitons in our study are formed as 
eigenstates of the nonlinear equation 10, one needs to 
demonstrate their stability under perturbation. First we 
suddenly modify the nonlinearity n of Eq. (jSJ after the 
dark soliton is obtained and study the resultant dynamics 
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for the dark soliton of Fig. 2 (a). In the case of the dark 
soliton of Fig. 2 (a) calculated using our scheme based on 
Eq. 1)1 0|t we suddenly jump the nonlinearity n from 100 to 
120 after the soliton is formed and observe its dynamics 
for 50000 units of time. The resultant dynamics is shown 
in Fig. 3. 

From Fig. 3 we find that even after giving a per- 
turbation by changing n, the resultant stationary dark 
soliton remains stable for a large interval of time (50000 
units of time) performing small breathing oscillation dur- 
ing which the central notch or minimum of the station- 
ary dark soliton remains absolutely stable at a; = 0. 
From Figs. 2 we find that in the iterative time evolu- 
tion method based on ansatz (JHJ the dark soliton may 
develop dynamical instability with the central notch ex- 
ecuting quasi-periodic oscillation around x = on time 
evolution before being destroyed without any perturba- 
tion whatsoever in a much smaller interval of time than 
that considered in Fig. 3. 

In addition to the perturbation studied in Fig. 3 we 
make two different types of perturbation to strengthen 
our claim of stability. First we study the dynamics by 
increasing the strength of the harmonic trap by 20%: (i) 
x 2 — > 1.2x 2 . These perturbations are symmetric around 
x = which may not displace the notch in the dark 
soliton from x = 0. We consider also the asymmetric 
perturbation (ii) u(x) — > u(x) + 0.02 x abs(u(x)) where 
abs denotes absolute value. As u(x) for the dark soliton 
is antisymmetric around x = 0, perturbation (ii) destroys 
the symmetry around x — 0. From Figs. 2 and 3 it is 
realized that it would be more difficult to have stability 
for a large nonlinearity. Hence in the next two studies 
on stability we consider only n = 100. The remarkable 
stability of the soliton under perturbations (i) and (ii) 
above is illustrated in Figs. 4 (a) and (b). Of these even 
the asymmetric perturbation (ii) leaves the central notch 
of the stationary dark soliton essentially stable at x — 0. 
The dynamics at small times in Fig. 4 (b) shows some 
asymmetry, which, instead of increasing, disappears at 
large times. From Fig. 2 (a) we find that the dark soliton 
calculated using ansatz © is destroyed after a small time 
of t = 100 without any perturbation whatsoever. 

There seems to be one difference between the behav- 
ior of the dark soliton in Figs. 4 and the corresponding 
behavior of a BEC in a trap. With a sudden change of 
the trap frequency, a BEC executes breathing oscillation 
where the central density fluctuates [2l]]. In Fig. 4 (a) 
we also observe a similar breathing oscillation. However, 
in this case because of the nature of the black soliton the 
central density remains zero, at least for a small change 
of trap frequency. In Fig. 4 (b) the dark soliton is given 
a asymetric displacement initially. Because of the sta- 
tionary nature of the initial dark soliton, the initial os- 
cillation tends to disappear after some time as the dark 
soliton settles to a stationary configuration. The scenario 
in Figs. 4 should change upon a large perturbation. 
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(a) „ (b) 

FIG. 4: The dark soliton |w(a;)| of Fig. 2 (a) at times 
t = 10000,20000,30000,40000,50000 calculated using the 
present scheme with input 1101 (a) as the harmonic oscillator 
potential x 2 is increased by 20% : x 2 — > 1.2a; 2 and (b) under 
the change u(x) —> u(x) + 0.02 x abs(u(a;)) at t = 0. 



III. HARMONIC PLUS OPTICAL-LATTICE 
TRAPS 

Now we consider the stationary dark soliton in a har- 
monic plus optical-lattice traps. A periodic optical- 
lattice trap is usually generated by a standing-wave laser 
beam of wave length A. In experiments the following su- 
perposition of a harmonic plus optical-lattice traps has 

been used [H El : 

V(x) = kx 2 + V sm 2 (2nx/X). (11) 

Here Vq is the strength of the optical-lattice potential. 
We have introduced a parameter k to control the strength 
of the harmonic potential. 

The search for a stationary dark soliton in potential 
(|ll|l is performed by introducing this potential in Eq. 
©. For this purpose we use the input u\(x,0) of Eq. 
(|10|l in Eq. JSJ with n = Vq = and perform time 
evolution using the Crank-Nicholson scheme [2l| . Again 
the discretization was performed with a space step 0.05 
and time step 0.0025 except for the calculation reported 
in Fig. 6 where we had to take a space step 0.02 and 
time step 0.0004. In the course of time evolution the 
appropriate nonlinearity and the optical-lattice potential 
are switched on slowly. Then the time evolution of the 
resultant equation is carried on until a converged solution 
is obtained. The results of the calculation for different n 
are shown in Figs. 5. These results are stationary and 
do not change with time evolution in the interval t = 
to t = 100000. 

However, it was more difficult to obtain convergence 
when nonlinearity n or strength Vq increases past 40. For 
small n(= 1, 10) convergence could be easily obtained for 
Vb up to 80 or so. For larger n(= 20) we could obtain 
convergence for Vb up to 40 or so. When both n and Vq 
are increased, smaller space and time steps are needed for 
convergence. This is understandable as a large nonlinear- 
ity with a strong optical-lattice potential could seriously 
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FIG. 5: The stationary dark soliton |w(a;)| of Eq. JSJ with 
potential (li lt with k = A = 1 vs. x for n = (a) —3 and 
(b) 10 with strength Vo of potential lilt given in respective 
figures. The plotted wave function remains stationary for time 
evolution of Eq. © over 100000 units of time t. 
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FIG. 6: The probability density |it(a:)| 2 of the stationary dark 
(black) soliton of Eq. @ with potential ill It with k — 0.001, 
Vo = 1.5 and A = 5 vs. i at different times. The dark soliton 
with a zero at x = remains stable for a long time time 
without executing quasi-periodic oscillation. 



jeopardize the numerical accuracy. There is no such dif- 
ficulty if the optical-lattice potential is removed. If A is 
reduced there is no convergence difficulty, so long as a 
finer discretization mesh is used. 

The possibility of generating dark solitons in a BEC 
with a harmonic plus optical-lattice traps l|llfl has also 
been investigated by other authors. Kevrekidis et al. 0] 
obtained stable solitons for their parameters (chemical 
potential, A etc.) for a very weak trap V(x) — kx 2 with 
k < 0.01. In this study, we could obtain such solutions 
for our parameter values including k = 1. 

In view of the present study, it seems that the insta- 
bility noted in Ref. could be due to the use of an 
initial non-stationary state, e. g. JHJ|, to generate the fi- 
nal stationary dark soliton. To substantiate our claim we 



consider a specific case highlighted by Kevrekidis et al. 
|13| as an example of instability of the dark soliton in an 
harmonic plus optical-lattice traps, e.g., for V(x) = kx 2 
with k = 0.001, V = 1.5, and A = 5 in Eq. (fTTfl . 

We repeated the above calculation of Kevrekidis et 
al. 0] using our approach. Because of the small value 
of k, the size of the condensate is much larger in this 
case compared to the condensates studied above and we 
had to take a much larger number of mesh points which 
results in large computing time and slower convergence. 
However, no quasi-periodic oscillation of the position of 
the minimum of the dark soliton was observed at large 
times. The present result is shown in Fig. 6 at differ- 
ent times. The central density in our calculation remains 
strictly zero over large time scales (t > 2000), whereas 
the calculation of Kevrekidis et al. ^| becomes un- 
stable for t > 500 with the notch of the dark soliton 
executing quasi-periodic oscillation. There is a difference 
in the two results, however, which prohibits us to make 
a quantitative comparison of the two calculations. In the 
present work the the wave function is normalized to unity, 
whereas in Ref. [13| the chemical potential has been fixed 
to unity. However, here we did not use a random pertur- 
bation to test the stability of the black soliton. 



IV. BLACK SOLITON IN THREE DIMENSIONS 

To further fortify the claim of numerical stability of 
our scheme we next test it in an axially-symmetric three- 
dimensional BEC under harmonic as well as optical- 
lattice traps. We consider the following GP equation for 
the BEC wave function ip(x, y; t) = <p(x, y, t)/x at radial 
position x, axial position y and time t: |25j 



.d d 2 1, a 2 2n_J_ 

dt dx 2 x dx dy 2 4 ^ x 2 



4:TT 2 K 



■ COS 



2iry 



(x,y;t) 



(x,y;t) = 0, (12) 



with normalization 



/OC f oc 

dy xdx\iP{x,y;t)\ 2 = 1, (13) 
-oo JO 

where (x 2 + v 2 y 2 )/4: is the axial harmonic trap and 
(4tt 2 k/X 2 ) cos 2 (27ry/A) the optical-lattice potential. 
Here length and time are expressed in units of l/\/2 = 
y/h/ (2mw) and respectively, with u> the radial trap 
frequency, m the atomic mass, v the axial parameter, 
and n — 8ir\f2Noa/l the scaled nonlinearity, where No is 
the number of atoms and a the scattering length. The 
numerical solution is calculated by the Crank-Nicholson 
discretization scheme where we used as in Ref. |2lJ a 
space step of 0.1 in both radial and axial directions and 
a time step of 0.001. The stationary dark soliton we are 
looking for is a nonlinear extension of the following solu- 





FIG. 7: Contour plot of the wave function \ip(x, y)\ of the 
black soliton of Eq. 1121 with harmonic trap alone for (a) 
v = 1, (b) v = 0.4, and (c) v = 0.1, and n = 100. 



FIG. 8: Contour plot of the black-soliton wave function 
\ip(x,y)\ of Eq. 11211 with with k = 4, A = 1 and (a) v = 1, 
(b) v = 0.4, and (c) v = 0.1, and n = 100. 



tion of Eq. ijT^)) with n = k = 0: 

/ v v 3/4 

#k,J/)=(^— J xyexp[-(x 2 + vy 2 )/A], (14) 

and can be found in a regular fashion by time evolution 
of Eq. ljT2|) with n = k = with Eq. (|T4"|) as the initial 
solution. This solution has a notch at y = in the axial y 
direction. In the course of time evolution the nonlinearity 
n and the optical-lattice strength n are slowly introduced 
until the final values of these parameters are attained. 

In this study we take the atoms to be 87 Rb and a final 
nonlinearity n = 100 together with u> = 2ir x 90 Hz so 
that ~ 0.8 /im. First we study the generation of 

a stationary dark soliton in the harmonic trap alone by 
setting k = in Eq. H12[) f° r three values of v — 1,0.4 
and 0.1. Contour plots of the wave function ^{x^y) of the 
stationary dark soliton in these three cases are exhibited 
in Figs. 7 (a), (b), and (c), respectively, where the central 
notch appears prominently and does not oscillate upon 
time evaluation. 

Next we study the stationary dark soliton in the above 
problem in the presence of an optical-lattice potential 
with k — 4 and A = 1 in addition to the above harmonic 
potential. The contour plots of the stationary dark soli- 
ton in this case are shown in Figs. 8. In addition to 
the more prominent notch at the center signaling a sta- 
tionary dark soliton, there are periodic lines along the 
axial direction due to the optical-lattice potential. This 
should be contrasted with the similar wave function in 



one dimension exhibited in Figs. 5 and 6, where there 
is also periodic modulation of the wave function due to 
the optical-lattice potential. Needless to say that these 
dark solitons in three-dimensions are stationary solutions 
of the GP equation and do not exhibit instability upon 
time evolution. 



V. CONCLUSION 

We show that the stationary dark (black) soliton of 
a trapped zero-temperature BEC is actually a nonlinear 
extension of the first vibrational excitation of the linear 
problem obtained by setting n = in Eq. ©. Based 
on this, we suggest a time-evolution calculational scheme 
starting from the linear problem with the harmonic po- 
tential alone while the nonlinearity and any additional 
potential (such as the optical-lattice potential) are slowly 
introduced during time evolution. This results in a stable 
numerical scheme for the stationary dark soliton as dur- 
ing time evolution the system always passes through suc- 
cessive stationary eigenstates of the nonlinear GP equa- 
tion until the desired stationary dark soliton is obtained 
also as a stationary eigenstate of the final GP equation. 
The present approach is equally applicable to both re- 
pulsive and attractive interactions in one and three di- 
mensions and eliminates the so called dynamical insta- 
bility of the dark soliton in the presence of a trap and 
yields a robust stationary dark soliton. To demonstrate 



7 



the robustness of the present numerical scheme, we il- 
lustrate the stable prolonged breathing oscillation of the 
dark soliton upon the application of a perturbation. We 
have presented our results in a dimcnsionless time unit. 
In a typical experimental situation this unit corresponds 
to about 1 ms. The stability of the dark soliton during 
few thousand units of time as considered in this paper 
is sufficient for this purpose. Time scales beyond a few 
seconds are completely unrealistic in this context, since 



the condensate only lives for about 15 seconds, with the 
soliton predicted to have a much shorter 'lifespan' due 
to thermal and quantum fluctuations even at very low 
temperatures. 
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